This is the analysis for an online survey to gather information about the acceptability of form factor and placements of light loggers. The survey was conducted using Prolific on three days in 2023, each collecting from a different sample (USA, UK, Worldwide).
This is the second part of the analysis, focusing on inference.
Import
We will start off by importing the data and setting factor levels.
Code
#Load Hmisc librarylibrary(Hmisc)library(tidyverse)library(gtsummary)library(gt)library(ggtext)library(rlang)library(patchwork)library(ordinal)library(rlang)library(glue)library(broom)library(ggridges)#Read Data in the long formatload(file ="data/cleaned/LightLoggerFormFacto_Long_CLEANED_DATA_2024-07-02.RData")
Inferential statistics
The base hypothesis is that the ratings depend on the Wearing Position, with Sex and Sample as control variables. Further, the effect of the Wearing Position might be dependent on Sex.
The inferential analysis fits cumulative link mixed models (CLMM) to the data. The models are fitted using the ordinal package. The models are fitted for each rating separately. The base model is fitted with the following formula:
Where \(E(logit(y))\) is the expected value of the rating on a logit scale, \(\alpha_0\) is the intercept, \(\beta_{i1}\) is the effect size of Wearing Position, \(\beta_{i2}\) is the effect size of Sex, \(\beta_{i3}\) is the interaction effect size of Wearing Position and Sex, \(\beta_{i4}\) is the effect size of Sample, and \(\alpha_j\) is the random intercept for each participant.
Models are compared through likelihood ratio tests. The following models are used for comparison:
A model without the interaction of Wearing Position and Sex (m2), a model without Sex (m3), a model without Sample (m4), a model with only Wearing Position (m5), and a null model, with only random effects (m0).
\(Rating \sim WearingPosition + Sex + Sample + (1|Id)\) (m2)
The following comparisons are made: 1. Interaction of Sex and Wearing Position (m1 vs. m2) 2. Inclusion of Sex (m2 vs. m3) 3. Inclusion of Sample (m1 vs. m4) 4. Inclusion of Wearing Position (m5 vs. m0)
The p-values are adjusted using the Benjamini & Hochberg (1995) method, also called the false discovery rate (fdr) with an n=4 for four comparisons.
Setup
Code
#convert some variables to factordata_long$sample_location <-as.factor(data_long$sample_location)data_long$record_id <-as.factor(data_long$record_id)#function generationInfe <-function(parameter, data) { parameter <-enexpr(parameter)#formula generation formulas <-inject(c(#full formulam1 =!!parameter ~ wearing_position*demographics_sex + sample_location + (1|record_id),#excluding interaction from m1m2 =!!parameter ~ wearing_position + demographics_sex + sample_location + (1|record_id),#further excluding sex from m2m3 =!!parameter ~ wearing_position + sample_location + (1|record_id),#excluding sample location from m1m4 =!!parameter ~ wearing_position*demographics_sex + (1|record_id),#excluding sample location from m3m5 =!!parameter ~ wearing_position + (1|record_id),#Null modelm0 =!!parameter ~1+ (1|record_id) ))#Model generation Infe_table <-tibble(name =names(formulas),formula = formulas,model =map(name, \(x) {switch(x,m1 = ,m2 = ,m3 = ,m4 =clmm(formulas[[x]], data = data, nAGQ =10,subset = demographics_sex !="Other"),m5 = ,m0 =clmm(formulas[[x]], data = data, nAGQ =10) ) } ) ) Infe_table}# Function to perform anova on model pairsanova_models <-function(pair) { one <- pair[[1]] two <- pair[[2]]inject(anova(Models$model[[!!one]], Models$model[[!!two]]))}#function to perform the model comparisonsModel_Comparisons <-function(Models) {# Define the model comparisonsmodel_comparisons <-list(c(1, 2),c(2, 3),c(1, 4),c(5, 6))tibble(desc =c("interaction of sex and wearing position","inclusion of sex","inclusion of sample location","inclusion of wearing position"),p.value =map(model_comparisons, \(x) anova_models(x) %>% .$`Pr(>Chisq)`%>% .[[2]] ) %>%unlist() %>%p.adjust(method ="fdr", n =4) %>% {case_when( . <0.001~"<0.001",TRUE~as.character(round(., 3)) )})}#probability of different answers#setting up a functionpred <-function(eta, theta, levels, cat =1:(length(theta)+1), inv.link = plogis) { Theta <-c(-Inf, theta, Inf) eta <-c(reference =0, eta) table <-sapply(cat, function(j)inv.link(Theta[j+1] - eta) -inv.link(Theta[j] - eta) ) colnames(table) <- levels table}data_significance_matrix <-function(Model, subset ="none") { Models <-tibble(reference =levels(data_long$wearing_position),data =map(reference, \(x) data_long %>%mutate(wearing_position =fct_relevel(wearing_position, x)) ),model =map(data, \(x) switch(subset,none =clmm(Model$formula, data = x, nAGQ =10),clmm(Model$formula, data = x, nAGQ =10,subset = demographics_sex !="Other") ) ),tidy =map(model, \(x) tidy(x) %>%filter(str_detect(term, "wearing_position")) %>%select(term, p.value) %>%mutate(term =str_remove(term, "wearing_position")) ) )Models %>%select(-data, -model) %>%unnest(tidy) %>%rowwise() %>%mutate(p.value =p.adjust(p.value, method ="fdr", n =7),different = p.value <=0.05) %>%ungroup() %>%rbind(tibble(reference =levels(data_long$wearing_position),term = reference,p.value =NA,different =FALSE))}#McFadden´s Pseudo R2, values between 0.2 bis 0.4 are an excellent model fitMcFadden <-function(model, nullmodel) { ll.null <- nullmodel$logLik ll.proposed <- model$logLik1- ll.proposed/ll.null}#Formula to plot the significance matrixSignificance_matrix <-function(sig_matrix, sex =FALSE){ sig_matrix %>%mutate(across(c(reference, term), \(x) fct_relevel(x, levels(data_long$wearing_position)))) %>%ggplot(aes(x=reference, y = term, fill = different)) +geom_tile(col ="grey30") +scale_fill_manual(values =c("white", "skyblue4"))+coord_fixed() +theme_void() +labs(title =label(data_long[parameter]),subtitle =case_when(!sex ~"<span style='color:skyblue4'>Significant differences</span> due to wearing position", sex ~"Significant differences due to <span style='color:skyblue4'>wearing position</span> and <span style='color:red'>sex</span>" ) )+theme(legend.position ="none",axis.text.x =element_text(angle =90, hjust =1),axis.text.y =element_text(hjust =1),plot.title =element_textbox_simple(hjust =0,padding =margin(0.5, 0, 0.5, 0, "cm")),plot.subtitle =element_textbox_simple(),plot.title.position ="plot", )}#table of probabilitiesprob_table <-function(beta, title){pred(beta, Model$Theta, levels =levels(data_long[[parameter]])) %>%as_tibble(rownames ="Wearing Position") %>%mutate(`Wearing Position`=str_remove(`Wearing Position`, "wearing_position"),`Wearing Position`=str_replace(`Wearing Position`, "reference", "Chest Pin") ) %>%gt() %>%fmt_percent() %>% gt::tab_caption(md(glue(title)))}
The only parameter which is supported by the data is Wearing Position. Sex, its interaction with Wearing Position, and Sample are not supported by the data. The resulting model m5 is used for the following analysis.
#if one wanted to print predictions for different subject percentilesstDev <-as.numeric(attr(VarCorr(Model)$record_id, "stddev"))^2# 5th percentile:prob_table(Model$beta *qnorm(0.05) * stDev, "Predicted probabilities for the **5th percentile** for **{parameter}**")
Predicted probabilities for the 5th percentile for appearance
Wearing Position
Extremely unattractive
Very unattractive
Unattractive
Neutral
Attractive
Very attractive
Extremely attractive
Chest Pin
0.95%
1.31%
6.69%
24.58%
51.44%
11.64%
3.40%
Wrist
0.30%
0.42%
2.28%
10.69%
50.30%
25.95%
10.05%
Necklace
0.07%
0.10%
0.57%
2.96%
26.44%
38.33%
31.52%
Sleeve collar
0.02%
0.02%
0.14%
0.73%
8.38%
24.71%
66.01%
Collar pin
0.00%
0.00%
0.03%
0.14%
1.71%
6.92%
91.20%
Neck loop
0.00%
0.00%
0.03%
0.15%
1.83%
7.37%
90.61%
Hat pin
0.00%
0.00%
0.01%
0.06%
0.75%
3.21%
95.97%
Glasses
0.00%
0.00%
0.00%
0.01%
0.17%
0.73%
99.09%
Code
#average:prob_table(Model$beta, "Predicted probabilities for the **average person** for **{parameter}**")
Predicted probabilities for the average person for appearance
Wearing Position
Extremely unattractive
Very unattractive
Unattractive
Neutral
Attractive
Very attractive
Extremely attractive
Chest Pin
0.95%
1.31%
6.69%
24.58%
51.44%
11.64%
3.40%
Wrist
1.49%
2.03%
9.94%
30.93%
45.55%
7.88%
2.17%
Necklace
2.59%
3.44%
15.42%
36.90%
35.66%
4.74%
1.25%
Sleeve collar
4.49%
5.71%
22.40%
38.68%
25.25%
2.76%
0.71%
Collar pin
8.38%
9.71%
30.38%
34.37%
15.35%
1.45%
0.37%
Neck loop
8.17%
9.51%
30.09%
34.66%
15.70%
1.49%
0.38%
Hat pin
11.29%
12.21%
33.18%
30.35%
11.65%
1.05%
0.26%
Glasses
18.87%
17.08%
34.55%
21.96%
6.82%
0.58%
0.14%
Code
# 95th percentile:prob_table(Model$beta *qnorm(0.95) * stDev, "Predicted probabilities for the **95th percentile** for **{parameter}**")
Predicted probabilities for the 95th percentile for appearance
Wearing Position
Extremely unattractive
Very unattractive
Unattractive
Neutral
Attractive
Very attractive
Extremely attractive
Chest Pin
0.95%
1.31%
6.69%
24.58%
51.44%
11.64%
3.40%
Wrist
2.95%
3.89%
16.97%
37.78%
33.14%
4.18%
1.09%
Necklace
11.13%
12.08%
33.07%
30.57%
11.82%
1.07%
0.27%
Sleeve collar
34.56%
21.48%
28.40%
12.09%
3.15%
0.26%
0.06%
Collar pin
73.80%
13.38%
9.48%
2.67%
0.61%
0.05%
0.01%
Neck loop
72.41%
13.96%
10.06%
2.86%
0.65%
0.05%
0.01%
Hat pin
86.61%
7.37%
4.54%
1.19%
0.27%
0.02%
0.01%
Glasses
96.72%
1.89%
1.06%
0.26%
0.06%
0.00%
0.00%
Model diagnostics
Code
#random effect distribution with random effectsci <- Model$ranef+qnorm(0.975)*sqrt(Model$condVar) %o%c(-1,1)ord.re <-order(Model$ranef)ci <- ci[order(Model$ranef),]plot(1:145, Model$ranef[ord.re], axes=FALSE, ylim=range(ci), xlab="Subject", ylab="Subject difference")axis(1, at=1:145, labels=ord.re)axis(2)for(i in1:145) segments(i, ci[i,1], i, ci[i,2])abline(h=0, lty=2)
Cumulative Link Mixed Model fitted with the adaptive Gauss-Hermite
quadrature approximation with 10 quadrature points
formula: rating ~ context + (1 | record_id)
data: data_long2
link threshold nobs logLik AIC niter max.grad cond.H
logit flexible 5800 -9957.53 19937.05 1209(8342) 8.62e-03 3.6e+02
Random effects:
Groups Name Variance Std.Dev.
record_id (Intercept) 1.451 1.205
Number of groups: record_id 145
Coefficients:
Estimate Std. Error z value Pr(>|z|)
contextpublic -0.45146 0.07412 -6.091 1.12e-09 ***
contextexercise -0.44588 0.07606 -5.862 4.57e-09 ***
contextwork -0.83751 0.07543 -11.103 < 2e-16 ***
contextsocial -0.91921 0.07517 -12.229 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Threshold coefficients:
Estimate Std. Error z value
Extremely unlikely|Very Unlikely -2.3870 0.1188 -20.100
Very Unlikely|Unlikely -1.5906 0.1163 -13.674
Unlikely|Neutral -0.6642 0.1149 -5.782
Neutral|Likely 0.1109 0.1145 0.968
Likely|Very likely 1.5633 0.1167 13.391
Very likely|Extremely likely 3.1163 0.1288 24.187
Code
summary(Model4)
Cumulative Link Mixed Model fitted with the adaptive Gauss-Hermite
quadrature approximation with 10 quadrature points
formula: rating ~ wearing_position + (1 | record_id)
data: data_long2
link threshold nobs logLik AIC niter max.grad cond.H
logit flexible 5800 -9401.64 18831.28 2076(16279) 9.58e-03 4.0e+02
Random effects:
Groups Name Variance Std.Dev.
record_id (Intercept) 1.851 1.36
Number of groups: record_id 145
Coefficients:
Estimate Std. Error z value Pr(>|z|)
wearing_positionWrist -0.33834 0.09403 -3.598 0.000321 ***
wearing_positionNecklace -1.01892 0.09474 -10.755 < 2e-16 ***
wearing_positionSleeve collar -1.39982 0.09718 -14.405 < 2e-16 ***
wearing_positionCollar pin -1.62095 0.09572 -16.935 < 2e-16 ***
wearing_positionNeck loop -1.99074 0.09701 -20.522 < 2e-16 ***
wearing_positionHat pin -2.59863 0.09921 -26.194 < 2e-16 ***
wearing_positionGlasses -2.62203 0.10032 -26.136 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Threshold coefficients:
Estimate Std. Error z value
Extremely unlikely|Very Unlikely -3.5293 0.1392 -25.352
Very Unlikely|Unlikely -2.6479 0.1362 -19.441
Unlikely|Neutral -1.6022 0.1338 -11.972
Neutral|Likely -0.7050 0.1324 -5.323
Likely|Very likely 0.9436 0.1327 7.109
Very likely|Extremely likely 2.6047 0.1431 18.196
Code
#get the predicted probabilities depending on context, save as a tablecontext_preds <-pred(Model3$beta, Model3$Theta, levels =levels(data_long2$rating)) %>%as_tibble(rownames ="Parameters") %>%mutate(`Parameters`=str_remove(`Parameters`, "wearing_position|context"),`Parameters`=str_replace(`Parameters`, "reference", "home") )context_preds_table <-context_preds %>%gt() %>%fmt_percent(decimals =0) %>% gt::tab_caption(md("Likelihood of wear ratings depending on context")) %>% gt::data_color(columns =-Parameters, method ="numeric", colors = scales::col_numeric(palette ="viridis", domain =c(0,0.30)) ) %>%cols_width(everything() ~px(100))
Warning: Since gt v0.9.0, the `colors` argument has been deprecated.
• Please use the `fn` argument instead.
This warning is displayed once every 8 hours.
The only parameter which is supported by the data is Wearing Position. Sex, its interaction with Wearing Position, and Sample are not supported by the data. The resulting model m5 is used for the following analysis.
#if one wanted to print predictions for different subject percentilesstDev <-as.numeric(attr(VarCorr(Model)$record_id, "stddev"))^2# 5th percentile:prob_table(Model$beta *qnorm(0.05) * stDev, "Predicted probabilities for the **5th percentile** for **{parameter}**")
Predicted probabilities for the 5th percentile for wearlikelihoodrating_public
Wearing Position
Extremely unlikely
Very Unlikely
Unlikely
Neutral
Likely
Very likely
Extremely likely
Chest Pin
1.79%
2.39%
8.79%
15.05%
45.90%
20.56%
5.51%
Wrist
0.52%
0.72%
2.88%
5.98%
34.91%
38.18%
16.80%
Necklace
0.12%
0.16%
0.67%
1.49%
12.95%
36.99%
47.63%
Sleeve collar
0.01%
0.01%
0.04%
0.09%
0.95%
5.22%
93.68%
Collar pin
0.00%
0.00%
0.02%
0.04%
0.42%
2.39%
97.12%
Neck loop
0.00%
0.00%
0.01%
0.02%
0.23%
1.35%
98.38%
Hat pin
0.00%
0.00%
0.00%
0.01%
0.08%
0.45%
99.46%
Glasses
0.00%
0.00%
0.00%
0.00%
0.03%
0.17%
99.80%
Code
#average:prob_table(Model$beta, "Predicted probabilities for the **average person** for **{parameter}**")
Predicted probabilities for the average person for wearlikelihoodrating_public
Wearing Position
Extremely unlikely
Very Unlikely
Unlikely
Neutral
Likely
Very likely
Extremely likely
Chest Pin
1.79%
2.39%
8.79%
15.05%
45.90%
20.56%
5.51%
Wrist
2.63%
3.44%
12.02%
18.49%
44.19%
15.44%
3.79%
Necklace
4.17%
5.26%
16.81%
21.92%
38.96%
10.49%
2.39%
Sleeve collar
9.52%
10.59%
26.13%
22.96%
25.04%
4.76%
1.00%
Collar pin
12.01%
12.60%
28.12%
21.72%
21.05%
3.73%
0.77%
Neck loop
14.13%
14.11%
29.11%
20.49%
18.40%
3.12%
0.64%
Hat pin
18.97%
16.93%
29.78%
17.65%
14.00%
2.22%
0.45%
Glasses
24.23%
19.11%
28.99%
14.89%
10.80%
1.64%
0.33%
Code
# 95th percentile:prob_table(Model$beta *qnorm(0.95) * stDev, "Predicted probabilities for the **95th percentile** for **{parameter}**")
Predicted probabilities for the 95th percentile for wearlikelihoodrating_public
Wearing Position
Extremely unlikely
Very Unlikely
Unlikely
Neutral
Likely
Very likely
Extremely likely
Chest Pin
1.79%
2.39%
8.79%
15.05%
45.90%
20.56%
5.51%
Wrist
5.94%
7.19%
20.93%
23.37%
33.33%
7.59%
1.66%
Necklace
22.16%
18.35%
29.43%
15.93%
11.92%
1.84%
0.37%
Sleeve collar
82.27%
9.46%
5.70%
1.57%
0.86%
0.12%
0.02%
Collar pin
91.35%
4.84%
2.66%
0.70%
0.38%
0.05%
0.01%
Neck loop
95.01%
2.84%
1.51%
0.39%
0.21%
0.03%
0.01%
Hat pin
98.31%
0.98%
0.50%
0.13%
0.07%
0.01%
0.00%
Glasses
99.36%
0.37%
0.19%
0.05%
0.03%
0.00%
0.00%
Model diagnostics
Code
#random effect distribution with random effectsci <- Model$ranef+qnorm(0.975)*sqrt(Model$condVar) %o%c(-1,1)ord.re <-order(Model$ranef)ci <- ci[order(Model$ranef),]plot(1:145, Model$ranef[ord.re], axes=FALSE, ylim=range(ci), xlab="Subject", ylab="Subject difference")axis(1, at=1:145, labels=ord.re)axis(2)for(i in1:145) segments(i, ci[i,1], i, ci[i,2])abline(h=0, lty=2)
The only parameter which is supported by the data is Wearing Position. Sex, its interaction with Wearing Position, and Sample are not supported by the data. The resulting model m5 is used for the following analysis.
#if one wanted to print predictions for different subject percentilesstDev <-as.numeric(attr(VarCorr(Model)$record_id, "stddev"))^2# 5th percentile:prob_table(Model$beta *qnorm(0.05) * stDev, "Predicted probabilities for the **5th percentile** for **{parameter}**")
Predicted probabilities for the 5th percentile for wearlikelihoodrating_work
Wearing Position
Extremely unlikely
Very Unlikely
Unlikely
Neutral
Likely
Very likely
Extremely likely
Chest Pin
2.31%
4.17%
11.02%
19.44%
41.41%
16.82%
4.82%
Wrist
0.43%
0.82%
2.48%
5.94%
30.15%
38.47%
21.70%
Necklace
0.05%
0.10%
0.33%
0.84%
6.36%
23.53%
68.78%
Sleeve collar
0.00%
0.00%
0.01%
0.02%
0.15%
0.78%
99.04%
Collar pin
0.00%
0.00%
0.01%
0.02%
0.17%
0.88%
98.92%
Neck loop
0.00%
0.00%
0.00%
0.01%
0.04%
0.23%
99.71%
Hat pin
0.00%
0.00%
0.00%
0.00%
0.00%
0.01%
99.99%
Glasses
0.00%
0.00%
0.00%
0.00%
0.01%
0.06%
99.92%
Code
#average:prob_table(Model$beta, "Predicted probabilities for the **average person** for **{parameter}**")
Predicted probabilities for the average person for wearlikelihoodrating_work
Wearing Position
Extremely unlikely
Very Unlikely
Unlikely
Neutral
Likely
Very likely
Extremely likely
Chest Pin
2.31%
4.17%
11.02%
19.44%
41.41%
16.82%
4.82%
Wrist
3.62%
6.30%
15.29%
23.00%
36.98%
11.72%
3.09%
Necklace
6.20%
10.03%
21.01%
24.86%
28.91%
7.21%
1.78%
Sleeve collar
15.88%
19.74%
27.26%
19.51%
14.27%
2.71%
0.63%
Collar pin
15.46%
19.43%
27.24%
19.79%
14.63%
2.80%
0.65%
Neck loop
20.81%
22.70%
26.71%
16.46%
10.89%
1.97%
0.45%
Hat pin
38.16%
26.23%
20.31%
9.16%
5.09%
0.85%
0.19%
Glasses
27.35%
25.11%
24.70%
13.16%
7.98%
1.39%
0.32%
Code
# 95th percentile:prob_table(Model$beta *qnorm(0.95) * stDev, "Predicted probabilities for the **95th percentile** for **{parameter}**")
Predicted probabilities for the 95th percentile for wearlikelihoodrating_work
Wearing Position
Extremely unlikely
Very Unlikely
Unlikely
Neutral
Likely
Very likely
Extremely likely
Chest Pin
2.31%
4.17%
11.02%
19.44%
41.41%
16.82%
4.82%
Wrist
11.45%
16.03%
26.23%
22.51%
18.98%
3.89%
0.92%
Necklace
50.68%
24.39%
15.14%
6.00%
3.15%
0.52%
0.12%
Sleeve collar
97.96%
1.33%
0.47%
0.15%
0.07%
0.01%
0.00%
Collar pin
97.72%
1.49%
0.53%
0.17%
0.08%
0.01%
0.00%
Neck loop
99.39%
0.40%
0.14%
0.04%
0.02%
0.00%
0.00%
Hat pin
99.97%
0.02%
0.01%
0.00%
0.00%
0.00%
0.00%
Glasses
99.83%
0.11%
0.04%
0.01%
0.01%
0.00%
0.00%
Model diagnostics
Code
#random effect distribution with random effectsci <- Model$ranef+qnorm(0.975)*sqrt(Model$condVar) %o%c(-1,1)ord.re <-order(Model$ranef)ci <- ci[order(Model$ranef),]plot(1:145, Model$ranef[ord.re], axes=FALSE, ylim=range(ci), xlab="Subject", ylab="Subject difference")axis(1, at=1:145, labels=ord.re)axis(2)for(i in1:145) segments(i, ci[i,1], i, ci[i,2])abline(h=0, lty=2)
The only parameter which is supported by the data is Wearing Position. Sex, its interaction with Wearing Position, and Sample are not supported by the data. The resulting model m5 is used for the following analysis.
#if one wanted to print predictions for different subject percentilesstDev <-as.numeric(attr(VarCorr(Model)$record_id, "stddev"))^2# 5th percentile:prob_table(Model$beta *qnorm(0.05) * stDev, "Predicted probabilities for the **5th percentile** for **{parameter}**")
Predicted probabilities for the 5th percentile for wearlikelihoodrating_home
Wearing Position
Extremely unlikely
Very Unlikely
Unlikely
Neutral
Likely
Very likely
Extremely likely
Chest Pin
1.12%
1.21%
4.98%
11.91%
42.81%
29.26%
8.72%
Wrist
0.02%
0.02%
0.08%
0.24%
2.08%
11.40%
86.16%
Necklace
0.03%
0.03%
0.14%
0.41%
3.48%
17.39%
78.51%
Sleeve collar
0.00%
0.00%
0.00%
0.00%
0.04%
0.26%
99.69%
Collar pin
0.00%
0.00%
0.00%
0.01%
0.05%
0.32%
99.63%
Neck loop
0.00%
0.00%
0.00%
0.00%
0.00%
0.03%
99.96%
Hat pin
0.00%
0.00%
0.00%
0.00%
0.00%
0.00%
100.00%
Glasses
0.00%
0.00%
0.00%
0.00%
0.00%
0.00%
100.00%
Code
#average:prob_table(Model$beta, "Predicted probabilities for the **average person** for **{parameter}**")
Predicted probabilities for the average person for wearlikelihoodrating_home
Wearing Position
Extremely unlikely
Very Unlikely
Unlikely
Neutral
Likely
Very likely
Extremely likely
Chest Pin
1.12%
1.21%
4.98%
11.91%
42.81%
29.26%
8.72%
Wrist
2.41%
2.55%
9.76%
19.52%
43.90%
17.68%
4.18%
Necklace
2.19%
2.33%
9.00%
18.52%
44.36%
19.01%
4.60%
Sleeve collar
4.92%
4.95%
16.70%
25.62%
36.04%
9.73%
2.04%
Collar pin
4.76%
4.80%
16.33%
25.42%
36.55%
10.04%
2.11%
Neck loop
7.17%
6.87%
21.02%
26.89%
29.84%
6.83%
1.38%
Hat pin
21.53%
15.19%
29.01%
19.53%
12.28%
2.06%
0.39%
Glasses
14.03%
11.63%
27.63%
24.19%
18.45%
3.40%
0.66%
Code
# 95th percentile:prob_table(Model$beta *qnorm(0.95) * stDev, "Predicted probabilities for the **95th percentile** for **{parameter}**")
Predicted probabilities for the 95th percentile for wearlikelihoodrating_home
Wearing Position
Extremely unlikely
Very Unlikely
Unlikely
Neutral
Likely
Very likely
Extremely likely
Chest Pin
1.12%
1.21%
4.98%
11.91%
42.81%
29.26%
8.72%
Wrist
42.35%
18.49%
22.86%
10.23%
5.13%
0.78%
0.15%
Necklace
30.13%
17.57%
27.39%
15.00%
8.33%
1.33%
0.25%
Sleeve collar
97.44%
1.34%
0.85%
0.25%
0.11%
0.02%
0.00%
Collar pin
96.92%
1.60%
1.03%
0.30%
0.13%
0.02%
0.00%
Neck loop
99.69%
0.16%
0.10%
0.03%
0.01%
0.00%
0.00%
Hat pin
100.00%
0.00%
0.00%
0.00%
0.00%
0.00%
0.00%
Glasses
99.99%
0.00%
0.00%
0.00%
0.00%
0.00%
0.00%
Model diagnostics
Code
#random effect distribution with random effectsci <- Model$ranef+qnorm(0.975)*sqrt(Model$condVar) %o%c(-1,1)ord.re <-order(Model$ranef)ci <- ci[order(Model$ranef),]plot(1:145, Model$ranef[ord.re], axes=FALSE, ylim=range(ci), xlab="Subject", ylab="Subject difference")axis(1, at=1:145, labels=ord.re)axis(2)for(i in1:145) segments(i, ci[i,1], i, ci[i,2])abline(h=0, lty=2)
The data supports Wearing Position, and its interaction with Sex. Sample is not supported by the data. The resulting model m4 is used for the following analysis.
Model comparisons for wearlikelihoodrating_exercise
desc
p.value
interaction of sex and wearing position
0.051
inclusion of sex
0.61
inclusion of sample location
0.506
inclusion of wearing position
<0.001
The only parameter which is supported by the data is Wearing Position. Sex, its interaction with Wearing Position, and Sample are not supported by the data. The resulting model m5 is used for the following analysis.
#if one wanted to print predictions for different subject percentilesstDev <-as.numeric(attr(VarCorr(Model)$record_id, "stddev"))^2# 5th percentile:prob_table(Model$beta *qnorm(0.05) * stDev, "Predicted probabilities for the **5th percentile** for **{parameter}**")
Predicted probabilities for the 5th percentile for wearlikelihoodrating_exercise
Wearing Position
Extremely unlikely
Very Unlikely
Unlikely
Neutral
Likely
Very likely
Extremely likely
Chest Pin
2.61%
3.34%
10.63%
16.48%
38.94%
21.46%
6.52%
Wrist
3.40%
4.26%
13.00%
18.62%
37.83%
17.83%
5.06%
Necklace
0.00%
0.01%
0.03%
0.05%
0.38%
2.12%
97.41%
Sleeve collar
1.78%
2.32%
7.74%
13.18%
38.45%
27.16%
9.36%
Collar pin
0.05%
0.07%
0.28%
0.59%
3.99%
17.65%
77.36%
Neck loop
0.00%
0.01%
0.02%
0.05%
0.33%
1.84%
97.75%
Hat pin
0.00%
0.00%
0.00%
0.01%
0.07%
0.38%
99.54%
Glasses
0.00%
0.00%
0.00%
0.00%
0.01%
0.06%
99.93%
Code
#average:prob_table(Model$beta, "Predicted probabilities for the **average person** for **{parameter}**")
Predicted probabilities for the average person for wearlikelihoodrating_exercise
Wearing Position
Extremely unlikely
Very Unlikely
Unlikely
Neutral
Likely
Very likely
Extremely likely
Chest Pin
2.61%
3.34%
10.63%
16.48%
38.94%
21.46%
6.52%
Wrist
2.41%
3.09%
9.95%
15.78%
39.05%
22.67%
7.06%
Necklace
16.20%
15.12%
27.56%
19.18%
16.82%
4.16%
0.96%
Sleeve collar
2.95%
3.74%
11.68%
17.49%
38.58%
19.76%
5.81%
Collar pin
8.35%
9.34%
22.60%
22.35%
27.08%
8.26%
2.01%
Neck loop
16.83%
15.48%
27.67%
18.85%
16.26%
3.99%
0.92%
Hat pin
25.05%
19.04%
27.14%
14.78%
10.95%
2.47%
0.56%
Glasses
38.04%
21.12%
22.82%
9.90%
6.46%
1.37%
0.30%
Code
# 95th percentile:prob_table(Model$beta *qnorm(0.95) * stDev, "Predicted probabilities for the **95th percentile** for **{parameter}**")
Predicted probabilities for the 95th percentile for wearlikelihoodrating_exercise
Wearing Position
Extremely unlikely
Very Unlikely
Unlikely
Neutral
Likely
Very likely
Extremely likely
Chest Pin
2.61%
3.34%
10.63%
16.48%
38.94%
21.46%
6.52%
Wrist
2.01%
2.60%
8.57%
14.21%
38.88%
25.35%
8.37%
Necklace
93.53%
3.62%
1.92%
0.55%
0.30%
0.06%
0.01%
Sleeve collar
3.82%
4.75%
14.17%
19.50%
36.96%
16.29%
4.50%
Collar pin
56.80%
18.82%
15.07%
5.34%
3.18%
0.64%
0.14%
Neck loop
94.35%
3.17%
1.67%
0.48%
0.26%
0.05%
0.01%
Hat pin
98.80%
0.68%
0.35%
0.10%
0.05%
0.01%
0.00%
Glasses
99.83%
0.10%
0.05%
0.01%
0.01%
0.00%
0.00%
Model diagnostics
Code
#random effect distribution with random effectsci <- Model$ranef+qnorm(0.975)*sqrt(Model$condVar) %o%c(-1,1)ord.re <-order(Model$ranef)ci <- ci[order(Model$ranef),]plot(1:145, Model$ranef[ord.re], axes=FALSE, ylim=range(ci), xlab="Subject", ylab="Subject difference")axis(1, at=1:145, labels=ord.re)axis(2)for(i in1:145) segments(i, ci[i,1], i, ci[i,2])abline(h=0, lty=2)
The data supports Wearing Position, and its interaction with Sex. Sample is not supported by the data. The resulting model m4 is used for the following analysis.
The only parameter which is supported by the data is Wearing Position. Sex, its interaction with Wearing Position, and Sample are not supported by the data. The resulting model m5 is used for the following analysis.
#if one wanted to print predictions for different subject percentilesstDev <-as.numeric(attr(VarCorr(Model)$record_id, "stddev"))^2# 5th percentile:prob_table(Model$beta *qnorm(0.05) * stDev, "Predicted probabilities for the **5th percentile** for **{parameter}**")
Predicted probabilities for the 5th percentile for restrict
Wearing Position
Extremely
Significantly
Quite a bit
Moderately
Somewhat
Slightly
Not at all
Chest Pin
0.18%
0.36%
0.80%
1.64%
5.20%
19.27%
72.54%
Wrist
0.01%
0.02%
0.05%
0.11%
0.38%
1.84%
97.58%
Necklace
0.00%
0.01%
0.02%
0.04%
0.12%
0.59%
99.23%
Sleeve collar
0.01%
0.02%
0.05%
0.11%
0.36%
1.74%
97.72%
Collar pin
0.00%
0.01%
0.02%
0.03%
0.12%
0.58%
99.24%
Neck loop
0.00%
0.00%
0.01%
0.02%
0.05%
0.25%
99.67%
Hat pin
0.00%
0.00%
0.00%
0.01%
0.03%
0.14%
99.82%
Glasses
0.00%
0.00%
0.00%
0.00%
0.01%
0.06%
99.93%
Code
#average:prob_table(Model$beta, "Predicted probabilities for the **average person** for **{parameter}**")
Predicted probabilities for the average person for restrict
Wearing Position
Extremely
Significantly
Quite a bit
Moderately
Somewhat
Slightly
Not at all
Chest Pin
0.18%
0.36%
0.80%
1.64%
5.20%
19.27%
72.54%
Wrist
0.71%
1.41%
3.02%
5.77%
15.26%
33.91%
39.94%
Necklace
1.26%
2.48%
5.12%
9.16%
20.87%
34.10%
27.02%
Sleeve collar
0.73%
1.45%
3.10%
5.91%
15.55%
34.05%
39.21%
Collar pin
1.27%
2.49%
5.15%
9.21%
20.95%
34.07%
26.86%
Neck loop
1.93%
3.71%
7.38%
12.27%
24.22%
31.13%
19.37%
Hat pin
2.59%
4.89%
9.35%
14.56%
25.60%
27.92%
15.09%
Glasses
4.08%
7.37%
13.01%
17.81%
25.69%
22.05%
9.99%
Code
# 95th percentile:prob_table(Model$beta *qnorm(0.95) * stDev, "Predicted probabilities for the **95th percentile** for **{parameter}**")
Predicted probabilities for the 95th percentile for restrict
Wearing Position
Extremely
Significantly
Quite a bit
Moderately
Somewhat
Slightly
Not at all
Chest Pin
0.18%
0.36%
0.80%
1.64%
5.20%
19.27%
72.54%
Wrist
2.65%
5.00%
9.54%
14.76%
25.68%
27.61%
14.76%
Necklace
7.98%
12.89%
18.90%
20.12%
21.33%
13.62%
5.16%
Sleeve collar
2.81%
5.28%
9.97%
15.21%
25.82%
26.89%
14.02%
Collar pin
8.10%
13.04%
19.02%
20.12%
21.19%
13.45%
5.09%
Neck loop
16.93%
21.34%
22.55%
17.01%
13.22%
6.69%
2.26%
Hat pin
26.96%
25.93%
20.87%
12.64%
8.44%
3.89%
1.26%
Glasses
48.33%
25.66%
13.69%
6.46%
3.75%
1.60%
0.50%
Model diagnostics
Code
#random effect distribution with random effectsci <- Model$ranef+qnorm(0.975)*sqrt(Model$condVar) %o%c(-1,1)ord.re <-order(Model$ranef)ci <- ci[order(Model$ranef),]plot(1:145, Model$ranef[ord.re], axes=FALSE, ylim=range(ci), xlab="Subject", ylab="Subject difference")axis(1, at=1:145, labels=ord.re)axis(2)for(i in1:145) segments(i, ci[i,1], i, ci[i,2])abline(h=0, lty=2)
a <-theme(plot.margin =margin(5, 10,5,5,"pt"))b <-theme(plot.margin =margin(5, 5,5,10,"pt"))(P1 + a + P2 + b) / (P3 + a + P4 + b) / (P5 + a + P6 + b) / (P7 + a + P8 + b) +plot_annotation(tag_levels ="A") &theme(plot.tag =element_text(size =20))